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Abstract 


Lithium-ion battery performance is strongly influenced by the ionic conductivity of the 
electrolyte, which depends on the speed at which Li ions migrate across the cell and relates 
to their solvation structure. The choice of solvent can greatly impact both solvation and dif- 
fusivity of Li ions. We use first principles molecular dynamics to examine the solvation and 
diffusion of Li ions in the bulk organic solvents ethylene carbonate (EC), ethyl methyl car¬ 
bonate (EMC), and a mixture of EC/EMC. We find that Ei ions are solvated by either car¬ 
bonyl or ether oxygen atoms of the solvents and sometimes by the PEg anion. Ei+ prefers 
a tetrahedrally-coordinated first solvation shell regardless of which species are involved, with 
the specific preferred solvation structure dependent on the organic solvent. In addition, we 
calculate Ei diffusion coefficients in each electrolyte, finding slightly larger diffusivities in the 
linear carbonate EMC compared to the cyclic carbonate EC. The magnitude of the diffusion 
coefficient correlates with the strength of Ei+ solvation. Corresponding analysis for the PEg 
anion shows greater diffusivity associated with a weakly-bound, poorly defined first solvation 
shell. These results may be used to aid in the design of new electrolytes to improve Ei-ion 
battery performance. 


1 Introduction 

There is a growing need to replaee gasoline and other fossil fuels with environmentally-friendly 
alternative energy sourees.^However, many of these alternative energy sources such as solar, wind, 
waves, and geothermal energy require advances in storage technology in order to become practical. 
Li-ion batteries are convenient, portable energy storage devices, which are currently used to power 
many handheld consumer devices and electric vehicles. Main components of a Li-ion battery 
include the anode, cathode, and electrolyte. Carbonaceous materials like graphite are often used 
for the anode due to the low cost of carbon and the ease with which Li intercalates into the material. 
Lithium transition metal oxides are frequently used for the cathode. The electrolyte provides the 
medium through which Li ions diffuse between the anode and cathode and usually consists of a 
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lithium salt dissolved in either an organie liquid, ionie liquid, or gel polymer.l^ 

The ehoiee of the eleetrolyte ean influenee the overall performanee of the Li-ion battery.l^Many 
existing batteries use organie liquids to dissolve the Li salt. Among the most eommon organie liq¬ 
uids used in eommereial batteries today are various earbonates, ineluding ethylene earbonate (EC), 
ethyl methyl earbonate (EMC), and mixtures of these. EiPEg is a frequently used salt in the elee¬ 
trolyte as it exhibits high ionie eonduetivity. A good organie solvent will be able to dissolve a high 
eoneentration of salt, resulting in a high dieleetrie eoefheient. A low solvent viseosity faeilitates 
ionie transport.^ Typieally, eyelie earbonates like EC have a high dieleetrie eonstant, but also have 
high viseosity, while linear earbonates like EMC have lower viseosity, but also a low dieleetrie 
eonstant. Moreover, some organie liquids like EC have a melting point above room temperature, 
so that they are not liquids over the entire operating temperature range. In order to resolve these 
issues and optimize the viseosity, dieleetrie eonstant, and melting point of the eleetrolyte for bat¬ 
tery performanee, linear and eyelie earbonates are often mixed.^ EC is a eommonly used liquid 
for mixed eleetrolytes sinee it is known to form a proteetive layer known as the solid eleetrolyte 
interphase (SEI)l^^on graphitie anodes, whieh prevents exeessive eleetrolyte deeomposition and 
promotes reversible interealation into and out of the anode. 

Experiments have been eondueted to examine the effeet of the eleetrolyte eomposition on 
Li transport and solvation using different speetroseopie teehniques sueh as ETIR,® Raman,^II^®l 
nuelear-magnetie resonanee (NMR) eieetrospray ionization mass speetroseopy,^ neutron 

seattering,^ and X-ray diffraetion.^ Some of these experiments attempted to determine the eoor- 
dination number around the Ei+ in using different Ei salts, solvents, and eoneentrationsRe- 
ported eoordination numbers in these works range from ~2 to 5. However, there has been relatively 
little experimental eharaeterization of the solvent moleeule orientation around Ei+. One sueh study 
was performed by Cazzanelli et al, who determined a eoordination number of ~2 in a mixture of 
EC and propylene earbonate (PC) at different eoneentrationsAt high eoneentrations, they eon- 
eluded that Ei+ was “sandwiehed” between two ring solvent moleeules. Reeent NMR experiments 
have also shown that there is a preferenee for EC to solvate Ei+ over DMC in mixed EC/DMC sys- 
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terns and that the carbonyl oxygen atoms are involved in solvation.l^ Separate NMR experiments 
on transport properties have also been carried out, in which experimental diffusion coefficients 
for Li+ in different electrolytes were found to range from 1 — 8 x 10^^ cm^/s at 30°C.^ Mixed 
EC/EMC systems have Li+ diffusion coefficient determined to be between 1.5 — 4.5 x 10^^ cm^/s, 
depending on salt concentration.^ 

Theoretical work has also been performed to understand Li transport and solvation in various 
carbonate electrolytes. To date, many of these simulations have been performed using classical 
force fields.122^^ Some of these studies indicate a coordination number of 4 where the Li+ inter¬ 
acts with carbonyl oxygen atoms of the carbonate .l22} 124 | 27 | 2 8l mixed carbonate systems such as 
EC/DMC, it was found that both EC and DMC participate in solvating Li+.l2^Eurthermore, there 
was a greater affinity for Li+ to dissociate from its counter-ion in cyclic carbonates relative to lin¬ 
ear carbonates.^ However, classical potentials are limited in their transferability and their ability 
to describe charge transfer effects. These limitations are not present in first-principles methods, 
which treat the electrons quantum mechanically. Static quantum calculations using cluster models 
have been used to study the energetics of different solvation structures,l^^^but these studies do 
not include the effect of the overall liquid environment. 

Only recently have first principles molecular dynamics (EPMD) based on forces from den¬ 
sity functional theory been used to study the solvation and transport properties of Li+ in differ¬ 
ent electrolytes Use of EPMD can be more predictive than classical force fields due to its 
parameter-free nature. EPMD has better transferability and can more accurately describe polariz¬ 
ability, charge transfer, and partial charges than classical potentials. It also provides a better basis 
for future comparative studies of the electrolyte solutions near interfaces and in reactive environ¬ 
ments. Previous work using EPMD by both Leung et a/.^^and Ganesh et a/.found Li+ solvation 
structures that agreed with previous classical force field simulations and were generally consistent 
with experiments. In addition, Ganesh calculated diffusion coefficients that were slightly higher 
than experimental values and previous theoretical work.l^ However, both of these works consid¬ 
ered systems with no more than ~300-400 atoms and simulation times of 13-25 ps, which may 
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have limited the ability to systematieally extraet quantities of interest. Both also used eertain ap¬ 
proximations, sueh as inereased hydrogen masses and large time steps, to make the simulations 
more eomputationally feasible. Therefore, it remains unelear how aeeurately these simulations 
fully deseribe Li+ solvation and diffusion in real systems. 

In this work, we earry out FPMD simulations to study solvation and diffusion in several ex¬ 
perimentally relevant earbonate-based organie Li-ion battery eleetrolytes. We eompare the eyelie 
earbonate solvent EC to the linear earbonate EMC, as well as an EC/EMC mixture. The typieal 
EiPEg salt is ehosen for this work. Solvation and transport properties for the EC/EiPp6 system have 
been previously studied with PPMD,l^2Hlllbut similar linear earbonate and eyelie/linear earbonate 
mixtures have been explored only using elassieal foree fields We analyze solvation struetures of 
Ei+ to examine how the ehoiee of solvent influenees the strueture and explore the range of possi¬ 
ble solvation struetures in eaeh solvent. Eurthermore, we examine the interaetion of the Ei+ and 
PEg and eompare their solvation properties. We ealeulate the diffusion eoeffieient in the different 
eleetrolytes to understand why Ei+ diffuses faster in one solvent than another and find eorrelations 
between solvation and diffusivity. The relation of solvation and diffusivity of PE6 is also studied 
and eompared to that of Ei+. Einally, in addition to performing these studies with larger system 
sizes and longer time seales than previous EPMD simulations, we further employ the ReaxFF foree 
field^^to quantify the effeets of finite size and time seales on the observed Ei+ solvation struetures 
and diffusivity. We antieipate that our findings ean be used to design new eleetrolytes that will 
improve the eyeling rate in batteries by tuning solvation to enhanee diffusivity. 

2 Computational Details 

We perform first prineiples moleeular dynamies using density funetional theory (DET) with the 
projeetor augmented wave (PAW) method^^^^^ and the PBE generalized gradient approximation 
exehange-eorrelation funetional,l ^^ l ^^ l as implemented in the VASF®^ software paekage. A 450 
eV plane-wave eutoff was used with Brillouin zone sampling restrieted to the E point. All molee- 
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ular dynamics simulations were performed in the NVT ensemble using a Nose-Hoover thermo- 

stat^SmH with the Nose frequeney of ~1000 em^^ eorresponding to a period of ~32 fs, and time 

step of 0.5 fs. Eaeh system was equilibrated for 5-7.5 ps at 330 K, followed by 30 ps of simulation 

time to gather statisties. A temperature of 330 K was used to mimie an intermediate Li-ion battery 

Fr^ 

operating temperature and to ensure that EC was not frozen = 310 K). 

Solvation struetures were eharaeterized with pair eorrelation funetions, ealeulated using a bin 
size of 0.03 A. Average eoordination numbers were eomputed from the integral of the pair eorre¬ 
lation funetion. We further quantify how tightly the ions are solvated by ealeulating the average 
residenee time of first shell solvent moleeules by fitting an exponential to the time eorrelation 
funetion 

^’soIyW = («(')■ H(0)>. (1) 

where H{t) is 1 if a given moleeule is within the first solvation shell and 0 otherwise.^^SHIllIlS 
distanee eutoff of the first solvation shell is taken from the first minimum in the pair eorrelation 
funetion between the solvated ion and the eenter of mass of eaeh solvent moleeule. A stretehed 
exponential of the form exp[—(t/r)^] gives the best fits, with x being the residenee time and /3 an 
adjustable parameter.l ^^ l ^^ l 

Diffusion eoeffieients for Li+ and PE^ were extraeted using two methods: (i) integration of the 
veloeity autoeorrelation funetion (VACE) via,^^^^^ 


D 


1 

3 



N 


Y^{vi{0)-vi{t))dt, 

i=l 


( 2 ) 


and (ii) analysis of the mean square displaeement (MSD) over time using the Stokes-Einstein 
relation,® 


6 At 


(3) 


Eor infinite statisties, using either the MSD or VACE to ealeulate the diffusion eoeffieient should 
produee the exaet same answer as they are mathematieally equivalent. However, for finite statisties. 
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they are not numerieally identieal,!^ so we ealeulate the diffusion eoeffieient using both methods. 
The VACF and MSD were both ealeulated by averaging over multiple trajeetory windows spanning 
the entire trajeetory with starting eonfigurations every 50 fs, using various window lengths from 5- 
15 ps in inerements of 2.5 ps. Equation|^was used to ealeulate D from the VACF for eaeh window 
length. Likewise, the slope of the linear regime in the MSD was used to ealeulate D from Eq.j^for 
eaeh window length. In eaeh ease, the values for D for eaeh window length were averaged to get 
a final estimate of the diffusion eoeffieient. Reported uneertainties refleet the standard deviation in 
this average. 

3 Results and Discussion 

3.1 Li+ Solvation in Ethylene Carbonate 

We performed first prineiples moleeular dynamies simulations of a single LiPFg moleeule dis¬ 
solved in a periodie box of 63 EC moleeules. This system eorresponds to 638 total atoms, with a 
density 1.32 g/ee and Li eoneentration of 0.23 M. Typieal Li eoneentration in eommereial batteries 
is 4-5 times larger, but this eoneentration was used to foeus on the solvation of a single Li+ with¬ 
out the effeets of other salt moleeules being present. We earried out two independent simulations 
where the Li+ and PF^ ions were initially either assoeiated or dissoeiated. In Fig. we show 
the Li-P distanee over eaeh of the two trajeetories, whieh indieate that the LiPFg remains either 
assoeiated or dissoeiated for the entire simulation. In addition, we display the trajeetories of the 
Li+ ion and the P from the PF^ ion where the eolor gradients (dark to light) indieate time. We 
find that Li+ and PF^ follow similar trajeetories even when dissoeiated, with separation ~5-8 A, 
giving evidenee of eorrelated motion. 

The pair eorrelation funetions between the Li ion and either the earbonyl oxygen atoms from 
EC, designated Oc, or the ether oxygen atoms, designated O^, are shown in Fig.|^ We find two dif¬ 
ferent solvation struetures, one where the PF^ stays apart from Li+ [Fig.j^a)] and another where 
they remain elose [Fig.|^b)]. For the ease where they remain assoeiated, we observe an average 
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Figure 1: (Left) Trajeetories for Li+ and P atom of PFg for 63 EC + 1 LiPF6 system. Color gradi¬ 
ents designate time, where darker eolors are earlier in time. (Right) Li-P distanee as a funetion of 
time for trajeetories that are initially dissoeiated or assoeiated. 

coordination number of ~2.5 for Li-Oc and ~0.5 for Li-0^ in the first solvation shell. In this 
case, the PFg also occupies a site in the first solvation shell with a total coordination number of 4. 
We label this solvation structure “3carbonylPF6.” In the case where Li and PFg are dissociated, we 
calculate a coordination number of 4 for Li-Oc- We denote this solvation structure “4carbonyl.” 
Representative snapshots of the “3carbonylPL6” and “4carbonyl” solvation structures are shown 
in Lig. The peak of the Li-Oc pair correlation function at approximately ~1.9 A agrees with 
previous classical^^^J^^and LPMdI^^^^ simulations which range between 1.7-2.0 A. Our total coor¬ 
dination number of 4 also agrees with previous theoretical world^^l ^ l ^^ l ^^ l ^^^ and experiments.^ 
We compare the thermodynamic stability of these two structures by computing the average relative 
energies over the trajectories, which are also indicated in Lig.[^ We find the “4carbonyl” solvation 
structure to be favorable by ~0.2 eV. We also examined the orientation of the solvent molecules 
around Li+ in each case by tracking the Oc-Li-Oc angle. Lig. [sj^a) shows the histogram of the 
Oc-Li-Oc angles for both the “3carbonylPL6” and “4carbonyr’ solvation structures. We see that 
both are peaked at ~110°, indicating a preference for a tetrahedral arrangement in both cases. We 
also find that the carbonyl group of the EC molecule tends to point toward the Li+, with Li-Oc-Cc 
angle ~140° as shown in Lig.j^b). Previous LPMD simulations by Ganeshl^have also determined 
this angle to be 140°, which is consistent with the experimental value of 138°,!^ but lower than 
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the classical force field valued of 150°. The tetrahedral pattern we observe for “4carbonyl” agrees 
with previously calculated solvation structures for Li+ in 

for “3carbonylPF6” suggests that regardless of the composition of the first solvation shell, Li+ 
prefers to be solvated in a tetrahedral fashion. 
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Figure 2: Li-Oc and Li-0^ pair correlation functions (solid lines) and their integrals (dashed 
lines) for (a) “4carbonyr’ and (b) “3carbonylPF6” solvation structures of EC. Snapshot of solvation 
structure and average energy relative to the lowest energy structure are shown in the inset. 



Figure 3: Histogram of a representative (a) Oc-Li-Oc angle and (b) Li-Oc-Cc angle (denoted in 
the insets) for the “3carbonylPF6” and “4carbonyr’ solvation structures of EC during the trajectory. 
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3.2 Li+ Solvation in Ethyl Methyl Carbonate 


Unlike EC, ethyl methyl earbonate (EMC) has a lower viseosity and freezes at a mueh lower 
temperature, but it also has a lower dieleetrie eonstant. We examined how the solvation struetures 
and ionie motion differ when a linear earbonate sueh as EMC is used as the solvent. The simulation 
system eonsisted of 42 EMC moleeules and 1 EiPE 5 (638 total atoms) with a density of 1.01 g/ee 
and eoneentration of 0.22 M. Again, we ran simulations with the EiPEg initially either assoeiated 
or dissoeiated. We found that the initially dissoeiated EiPEg re-assoeiated within 15 ps during 
the simulation, as shown in the Ei-P distanee plot in Pig. On the other hand, when EiPPg 
started assoeiated, it remained assoeiated throughout the simulation. Upon examination of their 
trajeetories, we found that when Ei+ and PP^ are initially dissoeiated, the PP^ moves toward and 
finds the Ei+. Overall, PP^ is observed to migrate further than Ei+. This is likely due to the Ei^ 
being more tightly solvated than the PPg , as diseussed later. 



Pigure 4: (Eeft) Trajeetories for Ei+ and P atom of PP^ for 42 EMC + 1 EiPPe system. Color gra¬ 
dients designate time, where darker eolors are earlier in time. (Right) Ei-P distanee as a funetion 
of time for trajeetories that are initially dissoeiated or assoeiated. 


We observe three different solvation struetures during the simulations. The pair eorrelation 
funetions of Ei-Oc and Ei-0^ for eaeh of these solvation struetures along with representative 
snapshots and relative energies are shown in Pig. We find that the lowest energy strueture is 
one where the Ei+ and PPg are assoeiated [Pig.j^a)]. This is a signiheant differenee eompared to 
EC. The total eoordination number of Ei^ in this preferred strueture is 4, eonsisting of 3 earbonyl 
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oxygens from 3 EMC molecules and the nearby PF^ . We designate this solvation structure as 
“3carbonylPF6." We note that a similar structure was also observed with EC, but it is not the lowest 
energy configuration. There is another solvation structure in EMC that is 4-fold coordinated, and 
this consists of 4 EMC molecules oriented such that Ei+ is solvated by 3 carbonyl oxygen atoms 
and one ether oxygen atom. In this case, the Ei+ and PF^ remain apart; we denote this structure as 
“3carbonylether" [Fig.|^b)]. The orientations of these two 4-fold coordinated solvation structures 
are analyzed in Fig. [^a), where we plot the Oc-Li-Oc bond angle for the “3carbonylPF6” and 
“3carbonylether” cases. For both these cases, there is a peak near 110°, similar to EC, indicating 
a preferred tetrahedral arrangement of the solvent molecules regardless of whether EMC or PF^ 
is solvating Ei+. We also see a preference for the carbonyl oxygen to point toward Ei+, from the 
Ei-Oc-Cc angle of ~150° shown in Fig. |^b). Similar conclusions were drawn for EC, but the 
Ei-Oc-Cc angle is slightly bigger for EMC than EC, which is also consistent with trends seen by 
Borodin et al.^ using classical force fields. Unlike for EC, however, for EMC we also found a 
non-tetrahedral structure where the Ei+ is solvated by 4 total ether oxygen atoms belonging to two 
EMC molecules (2 ether oxygen atoms per EMC) and PFg , in a square pyramidal-like fashion. 
We label this structure “4etherPF6” [Fig.|^c)]. This structure is not energetically preferred, being 
more than 0.4 eV higher in energy than “3carbonylPF6.” For EMC, a corresponding “4carbonyl” 
solvation structure, which is favored by EC, is not observed in either of the trajectories. Steric 
issues likely prohibit this structure for EMC, since the length of the molecule makes it unfavorable 
to have four EMC molecules around Ei+. With PF^ similar in size to EC, the “3carbonylPF6” 
structure is preferred for EMC, instead of the “4carbonyl” structure as for EC. 

3.3 Li+ Solvation in 3:7 Mixture of EC/EMC 

We also examined the effect of mixing different organic solvents for the electrolyte. We studied a 
mixture of EC and EMC in a 3:7 ratio, mimicking previous experiments.^^ The simulation system 
consisted of 15 EC molecules, 35 EMC molecules, and 1 EiPFg (683 total atoms) with a density 
of 1.165 g/cc and a concentration of 0.23 M. As in the previous studies, we began the simulations 


11 


EMC “3carbonylPF6” EMC “3carbonylether” EMC “4etherPF6” 



Figure 5: Li-Oc and Li-0^ pair correlation functions (solid lines) and their integrals (dashed 
lines) for (a) “3carbonylPF6,” (b) “Scarbonylether,” and (c) “4etherPF6” solvation structure of 
EMC. Snapshot of solvation structure and average energy relative to the lowest energy structure 
are shown in the inset. 



Figure 6: Histogram of a representative (a) Oc-Li-Oc angle and (b) Li-Oc-Cc angle (denoted 
in the insets) for “3carbonylPF6” and “3carbonylether” solvation structures of EMC during the 
trajectory. 
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with the LiPF 5 either assoeiated or dissoeiated. In Fig. |7| we plot the Li-P distanee for both eases 
and show the trajeetories of the two ions. We see that the ions remain assoeiated or dissoeiated for 
the duration of the simulation. In addition, based on the trajeetories, we see that the Li and P in 
both the assoeiated and dissoeiated eases follow very similar paths, indieative that the ions behave 
more like in EC than in EMC, with a large degree of eorrelated motion even when the ions are 
separated. This similarity to EC oeeurs even with the mixture eontaining only 30% EC. Therefore, 
adding only a small pereentage of EC to the system results in a dramatie ehange in the ion motion 
with respeet to one another. 



Eigure 7: (Left) Trajeetories for Li+ and P atom of PEg for 15 EC + 35 EMC + 1 LiPEg system. 
Color gradients designate time, where darker eolors are earlier in time. (Right) Li-P distanee as a 
funetion of time for trajeetories that are initially dissoeiated or assoeiated. 

Prom these trajeetories, we find three different solvation struetures: two struetures where the 
Li+ and PP^ ions stay apart and one strueture where they remain together. Snapshots of eaeh sol¬ 
vation strueture, their pair eorrelation funetions for Li-Oc and Li-0^:, and their relative energies 
are shown in Pig. Por the two solvation struetures with LiPPg dissoeiated, one shows a eoordi- 
nation number of 4 for Li-Oc, indieating that 4 earbonyl oxygen atoms solvate Li+; the other has 
a eoordination number of ~3 for Li-Oc and ~1 for Li-Oc, where 3 earbonyl oxygen atoms and 1 
ether oxygen atom solvate Li+. We denote these solvation struetures “4earbonyl” [Pig.j^a)] and 
“3earbonylether” [Pig.j^b)], respeetively. Por the one strueture where LiPPg is assoeiated, there is 
a eoordination number of 2 for Li-Oc and ~1 for Li-Oc, indieating 2 earbonyl oxygen atoms and 
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1 ether oxygen atom solvate Li+ along with the PF^. We label this strueture “2earbonyletherPF6” 
[Fig. [^e)]. The “4earbonyl” and “Searbonylether” solvation struetures are nearly energetieally 
equivalent with only a 0.04 eV differenee between the two, whereas the “2earbonyletherPF6” is 
more than 0.3 eV higher in energy. In this 3:7 EC/EMC mixture, Li+ prefers to be separated from 
PFg , similar to pure EC. In addition, we observe that one EMC in the “3earbonylether” solvation 
structure is replaced by an EC molecule during one of the trajectories, leaving 2 EC molecules 
and 2 EMC molecules involved in the first solvation shell and forming the “4carbonyl” solvation 
structure. This result indicates a strong preference for EC molecules to solvate Ei+, considering 
the small fraction of EC in the system. 

Recent NMR experiments performed on mixed EC/DMC systems have also shown a strong 
preference for EC to solvate Ei+ as opposed to a linear carbonate such as DMC (or EMCj.l^This 
is consistent with our observation that one EMC molecule is replaced by an EC molecule in the first 
solvation during the course of our simulation. In addition, classical simulations also observe an 
equal amount of cyclic carbonate, EC, and linear carbonate, DMC, when the ions are dissociated.^ 
Therefore, our results are consistent with experiments and classical simulations and provide strong 
evidence that Ei+ prefers to be solvated by EC when present in these systems. As described above, 
this tendency likely is related to steric effects and the fact that EC and PE 5 have similar sizes. We 
also examined the orientation of the solvent molecules for all three cases. In Eig. |^a), we show 
histograms of the Oc-Ei-Oc angle for each solvation structure. Again, we find the peak of the 
histogram for all three cases near 110°, with Ei+ preferring a tetrahedral solvation structure. In 
Eig. I^b), we plot the Ei-Oc-Cc angle, and find that the carbonyl oxygen atoms prefer to point 
toward Ei+ with an angle of ~140°, similar to the previous cases (and particularly similar to pure 
EC). 

Overall, our results show that even small variations in the organic solvent can dramatically 
change the preferred solvation structure, although Ei always prefers to be coordinated tetrahedrally 
in its first solvation shell regardless of which species are around it. 
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3:7 EC/EMC “4carbonyl” 3:7 EC/EMC “3carbonylether” 3:7 EC/EMC “2carbonyletherPF6” 



Figure 8: Li-Oc and Li-0^ pair correlation functions (solid lines) and their integrals (dashed 
lines) for (a) “4carbonyl,” (b) “Scarbonylether,” and (c) “2carbonyletherPF6” solvation structures 
of mixed EC/EMC. Snapshot of solvation structure and average energy relative to the lowest energy 
structure are shown. 


0.05 


^ 0.04 

*(73 

Q 0.03 


2 0.02 

Xi 

0.01 


2carbonyletherPF6 
3carbonylether 
4carbonyl 


a) 


3:7 EC/EMC 









40 60 80 
Oj.-Li-0^ Ang 


0.05 


— 2carbonyletherPF6 
^0.041- — 3carbonylether 

— 4carbonyl 


3:7 EC/EMC 


C/3 

Q 0.03 


3 0.02 

I 0.011- 





b) 


20 140 160 180 
e (Degrees) 


JjJ 


40 60 80 100 12( 


M 


160 180 


Li-O^-C^ Angle (Degrees) 


Eigure 9: Histogram of a representative (a) Oc-Ei-Oc angle and (b) Ei-Oc-Cc angle (denoted in 
the insets) for the “4carbonyl,” “2carbonyletherPE6,” and “3carbonylether” solvation structures of 
mixed EC/EMC. 
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3.4 Solvation Structure of PF^ 

Thus far, we have primarily examined the solvation of the Li ion. We now turn our attention 
to the solvation strueture of the PF^ anion and compare it to that of Li+. There is relatively 
little discussion about the solvation structure of the counter-ion in previous FPMD and classical 
simulations in the dissociative limit. The EC system is used here for illustration, although EMC 
and the mixture show similar phenomena. In Fig. [T^ we show the pair correlation functions of the 
center of mass of the EC molecules and either Li [Fig.[T0]^a)] or P [Fig.fT^b)], for both trajectories 
where the ions were initially either associated or dissociated. The Li-EC pair correlation functions 
are well structured with a sharp first-shell peak, indicating that Li+ has a well-defined solvation 
structure with a clear coordination number ~4 in the first solvation shell, as already discussed 
in detail. Conversely, the P-EC pair correlation functions are very broad, suggesting that many 
EC (solvent) molecules dynamically rotate in and out of the first solvation shell on a short time 
scale. To quantify the strength of solvation, we computed the average residence time of the first- 
shell solvent molecules around each of the ions. Eor PE^, the residence time of EC was 43- 
90 ps and of EMC was 24-29 ps. (The ranges correspond to variations for trajectories with LiPEg 
either associated or dissociated.) Eor Li+, the residence times were well beyond the length of the 
simulations for both EC and EMC, with the first solvation shells showing little solvent exchange 
during the trajectories. While previous EPMD simulations^^^^ for EC report no solvent exchange 
during their trajectories, we see occasional solvent exchange that preserves the tetrahedral solvation 
structure. Thus, we determine that PE^ is much more weakly solvated, with a poorly-structured 
solvation shell and very short solvent molecule residence times, than Li+. Also, we find that EC 
solvates PE^ somewhat better than EMC, based on the relative residence times. 

We have already shown in Eigs.[^|^ andthat PE^ appears to move farther than Li+, despite 
the much heavier mass of the anion. The increased diffusivity of PEg compared to Li+ is connected 
to the respective solvation structures and is discussed in detail in the following section. 
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Figure 10: (a) Li-EC and (b) P-EC pair correlation functions (solid lines) and their integrals 
(dashed lines) for EiPE 5 in EC, with the salt molecule initially either associated or dissociated. 

3.5 Li+ and PF^ Transport Properties 

The transport properties of Ei+ and PE^ in each of the electrolytes were studied by comparing 
their diffusion coefficients. Table [I] compares the Ei+ diffusion coefficients for each electrolyte, 
calculated using both the MSD and VACE methods (see Sec.j^for details). We focus on the cases 
with the solvation in the preferred configurations, which generally corresponds to EiPE6 being 
dissociated, as discussed above. Trends in the diffusion coefficients are generally consistent across 
the different systems regardless if MSD or VACE is used for the calculation. Statistical errors of 
the calculated values are on the order of 1 — 2 x 10^^ cm^/s as shown in Table [I] or as much as 
50%. Improving these errors is nontrivial as it would require running many (>20) independent 
simulations with one Ei ion or a single simulation with many more Ei ions (and correspondingly 
larger system size). 

Values for EC are approximately a factor of two different than previous EPMD results^^^^ 
which calculated the diffusion coefficient to be ~ 1 x 10 ^ cm^/s at temperatures between 310- 
450K, likely due to the large uncertainty associated with the short DET runs. The highest Ei 
diffusion is seen in EMC. This is in agreement with experiments^ and classical simulations^ that 
find faster diffusion in linear carbonates compared to cyclic carbonates. Surprisingly, the mixed 
electrolyte shows slower Ei+ diffusion than pure EC, whereas it would be expected to fall between 
EC and EMC. We suspect the reason why the diffusion coefficient of the mixture does not fall 
between EC and EMC is related to statistical error from the rather short DET simulation and the 
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variation of solvation structure during the trajeetory. However, we note that the similar solvation 
struetures of Li+ in EC and the EC/EMC mixture result in similar diffusivities in these eases, 
whieh are distinet from that in pure EMC. In faet, the first solvation shell of Ei+ in the EC/EMC 
mixture ean eontain up to 2 EC moleeules, as diseussed above so it is expeeted that the value of 
the diffusion eoeffieient would be eloser to that of pure EC. 

Table 1: Calculated Li+ diffusion coefficients in each electrolyte, from the slope of the mean- 
square displacement (MSB) and integral of the velocity autocorrelation function (VACF). For 
each electrolyte, the most stable solvation configuration(s) were considered. 



Ei+ Diffusion Coefficient 


(10- 

cm^/s) 

Electrolyte eomposition 

MSD 

VACF 

63 EC -1-1 EiPFe 

5.2 ±0.8 

7.9 ± 1.3 

42 EMC -1-1 EiPFe 

9.6 ± 1.6 

10.1 ± 2.1 

15 EC -1- 35 EMC -i- 1 EiPFe 

2.6 ± 1.3 

5.1 ± 1.1 


Furthermore, analysis of the range of trajeetories shows that slower diffusion tends to oeeur 
for oases where the Ei^ solvation struoture is more energetioally preferred. A similar observation 
is noted when the coordination number is greater than 4 \i.e., the non-tetrahedral EMC ease in 
Fig.|^o)]. These results reveal that the magnitude of the diffusion eoeffieient is strongly dependent 
on how tightly solvated the Ei+ is by its solvent moleeules. When the Ei+ is more tightly solvated, 
the diffusion eoeffieient is smaller than when it is weakly solvated. We oonolude that EC solvates 
Ei+ better than EMC, as indioated by the lower diffusion eoeffieient. 

Table shows the diffusion coeffieients for PF^ , ealeulated by traeking the P atom. Overall, 
the values are larger than for Ei+, consistent with the weaker solvation structure diseussed above. 
We also note the higher diffusivity in EMC eompared to EC, whieh is due to the weaker solvation 
by EMC as evideneed by the shorter first-shell solvent moleeule residenee time. 


3.6 Finite Size and Time Scale Effects 

Moleeular dynamies simulations based on traditional Kohn-Sham density funetional theory imple¬ 
mentations are limited to moderate system sizes on the order of hundreds of atoms and time scales 
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Table 2: Calculated PF^ diffusion coefficients in each electrolyte, from the slope of the mean- 
square displacement (MSD) and integral of the velocity autocorrelation function (VACF). For 
each electrolyte, the most stable solvation configuration(s) were considered. 



PF^ Diffusion Coefficient 


(10- 

^ cm^/s) 

Electrolyte composition 

MSD 

VACF 

63 EC -1- 1 EiPF6 

7.1 ±0.9 

9.2 ± 1.0 

42 EMC -1- 1 EiPF6 

30.8 ± 8.8 

28.6 ±5.7 

15 EC -1- 35 EMC -i- 1 EiPF6 

5.7 ± 2.4 

9.5 ± 1.4 


of 10s of picoseconds. In order to gauge finite size and time scale effects on the solvation struc¬ 
tures and diffusion coefficients that we calculated using DFT, we used the ReaxFF reactive force 
fieidl 3^F0| II] as implemented in the LAMMPS^^^^^ software package to run much longer molecular 
dynamics trajectories of 1 ns with system sizes up to ~6400 atoms. 

First, we assessed the quality of the ReaxFF force field to reproduce the results of DFT, using 
the EC electrolyte as a test case. These simulations were performed under NVT conditions using 
a Nose-Hoover chain thermostat with 3 Nose-Hoover chains and time step of 0.25 fs. The Nose 
frequency was set to ~1333 cm^^ corresponding to a period of ~25 fs. A system of 630 EC and 
10 EiPE 5 molecules was equilibrated for 125 ps at 330 K followed by 1 ns of simulation time. This 
test system contained the same concentration of EiPEg as the DPT simulation, but with 10 times 
more Ei ions to gather better statistics. 

In Pig. [TT]^a), we show the Ei-0 pair correlation function and its integral for both the ReaxFF 
trajectory and the DFT trajectory. We find that both methods predict the same coordination number 
of ~4 carbonyl oxygen atoms for the first solvation shell. Also, the Oc-Ei-Oc angle distributions 
are nearly identical, with a peak at ~110° indicating a tetrahedral arrangement, as shown for 
example in Fig. However, the second peak in the Ei-0 pair correlation function (associated 
with the ether oxygen atoms) is slightly different, with DFT predicting a broader peak centered 
further away than ReaxFF. This second peak is much sharper and more structured with ReaxFF. 
We can understand this difference more deeply by examining the Ei-Oc-Cc angle. Figure [TT]^b) 
shows that DFT exhibits a broad distribution of angles centered at ~140°, while ReaxFF predicts a 
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narrow distribution of angles around ~90°. Representative snapshots of these solvation struetures 
from ReaxFF and DFT are displayed in the insets of Fig.[TT]^b), whieh illustrate how the 90° angle 
from ReaxFF results in the ether oxygen atoms being eloser to Li+ than in the DFT simulations. 
In addition, the seeond Li-0 peak in the pair eorrelation funetion from ReaxFF is mueh sharper 
than from DFT, beeause this “bent” solvation strueture is more rigid, presumably from additional 
interaetions between Li and O^. In the DFT simulations, the earbonyl dipoles of the EC moleeules 
point toward the Li+, but also exhibit more rotational fluetuations, giving rise to the broadening of 
the seeond peak in the Li-0 pair eorrelation funetion and of the Li-Oc-Cc angular distribution. 

ReaxFF also prediets oeeasional EC dimerization over the eourse of the trajeetory where the 
earbonyl earbon of one EC moleeule internets with the earbonyl oxygen of another EC moleeule. 
This dimerization is ineonsistent with DET and the ehemieal inertness of the EC liquid, but only 
oeeurs for less than 5% of the moleeules over 1 ns. We further eompared the diffusion eoeffieients 
for both Ei+ and PE^ and found that the ReaxFF values are within ~40-50% of the DFT values. 
Thus, we do note differenees between DFT and ReaxFF for these systems, but the applieability of 
ReaxFF to study finite size and time effeets appears valid. 



r (Angstrom) (Degrees) 

Figure 11: (a) Ei-0 pair eorrelation funetions (solid lines) and their integrals (dashed lines) for 
dissoeiated EiPFg in EC, eomparing DFT and ReaxFF trajeetories. (b) Histogram indieating the 
distribution of the Ei-Oc-Cc angle for ReaxFF and DFT trajeetories. Insets show representative 
snapshots of the ReaxFF and DFT Ei+ solvation struetures. 

To determine the effeet of time seale on solvation strueture, we ealeulated the Ei-0 pair eorre- 
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lation function and its integral using a small 30 ps segment of the 630 EC + 10 LiPFg 1 ns ReaxFF 
trajectory. In Fig. [T^a), these results are compared to the pair correlation function obtained when 
the entire 1 ns trajectory is used. We see that peak locations, intensities, and coordination num¬ 
ber are very similar, indicating that the time scale used for the MD simulations does not have a 
significant effect on the solvation structure. Similarly, to determine the effect of finite size on 
the solvation structure, we ran a 63 EC -i- 1 FiPFg MD simulation with ReaxFF for 1 ns under 
NVT conditions at 330K. The Fi-0 pair correlation function and integral for this 63 EC -i-1 EiPFg 
system are compared to that calculated for the 630 EC -i- 10 EiPF6 system in Fig. [T^b). We see 
that although the peak intensities vary slightly, the peak locations and the integral of the pair cor¬ 
relation function are comparable. These results were confirmed with several additional system 
sizes as well. Based on these observations, we don’t expect the solvation structures calculated 
using DFT to change significantly when going to longer time scales or larger system sizes at fixed 
concentrations. 



Figure 12: (a) Comparison of Ei-0 pair correlation functions and their integrals for 30 ps trajectory 
and 1 ns trajectory, (b) Comparison of Ei-0 pair correlation functions and their integrals for a 63 
EC -I- 1 EiPFg system size and 630 EC -i- 10 EiPF6 system size, where the concentration of the two 
systems are fixed. 

We also examined the effects of time scale and finite size on the diffusion coefficients for both 
Ei+ and PF^ in pure EC using ReaxFF. We compared the diffusion coefficient for Ei+ and PFg , 
calculated using the slope of the mean square displacement, from a 30 ps segment of the 630 EC 
-I- 10 EiPF6 trajectory to that from the entire 1 ns trajectory. We found that using only 30 ps of 
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the trajectory resulted in a difference of ~46% for Li+ and 38% for PF^ as opposed to using the 
whole 1 ns trajectory in the calculation. Based on these results, we expect that our DFT-calculated 
diffusion coefficients may vary up to 50% by running longer simulations. Classical simulations 
have also been used previously to study size effects, which found a 10% difference in the diffusion 
coefficient between large (480 solvent molecules) and small (240 solvent molecules) box sizes 
However, more Li atoms were included in these simulations, which increased the statistical sam¬ 
pling and resulted in correspondingly smaller finite-size effects. We further compared the diffusion 
coefficients calculated from 1 ns trajectories of systems containing either 63 EC -i- 1 LiPF6 or 630 
EC -I- 10 LiPF6, to judge the effects of finite size without changing the statistical sampling. The 
diffusion coefficients between the smaller and larger systems differed by ~6% for Li+, but up to 
~48% for PFg . Although the discrepancy for PF^ is fairly large, it is comparable with the ap¬ 
proximately 50% uncertainty we find for the time scale effects. Tests on EMC and the EC/EMC 
mixture showed similar results as EC. 

To summarize, we find that the Ei+ solvation structures do not change significantly when as¬ 
sessing finite size and time scale effects, although there are some differences between the ReaxFF 
and DET solvation structures. Eurthermore, uncertainties up to ~50% in the diffusion coefficients 
are expected for simulations run at short time scales and smaller length scales. This uncertainty is 
still small enough to allow us to draw qualitative and semi-quantitative conclusions as above, given 
the relative statistical errors in our computed values. However, it will be important in the future to 
use DET at still larger length and time scales to further reduce uncertainties in solvation structures 
and diffusion coefficients. 


4 Conclusions 

We found multiple possible solvation structures of Ei+ in each of the electrolytes studied here, 
including ethyl carbonate (EC), ethyl methyl carbonate (EMC), and EC/EMC mixture. While 
previous literature on EC and EiPEg has focused on solvation by 4 carbonyl oxygen atoms, we 


22 


observed solvation struetures that also inelude PFg in the first solvation shell. We found that 
the preferred solvation strueture is strongly dependent on the solvent. For EC, Li+ prefers to 
be solvated by 4 earbonyl oxygen atoms and to remain dissoeiated from PFg , while for EMC it 
shows some preferenee for only 3 earbonyl oxygen atoms and to remain elose to the PE^ . The 
3:7 EC/EMC mixture shows a preferenee for Ei+ to remain dissoeiated from PE^ , but has two 
energetieally similar solvation struetures where Ei+ is solvated either by 4 earbonyl oxygen atoms 
or by 3 earbonyl oxygen atoms and 1 ether oxygen atom. In all oases, Ei+ prefers to be solvated 
in a tetrahedral arrangement, though this does not rule out the possibility of a slightly higher- 
energy non-tetrahedral strueture forming, as seen in the EMC “4etherPE6” ease. Comparisons of 
solvation struetures for Ei+ and PE^ reveal that Ei+ is more strongly solvated, assooiated with 
lower mobility in the eleetrolyte. Caloulations of first shell solvent moleoule residenoe times show 
that there is a slight preferenee for PE^ to be solvated by EC over EMC, although both show weak 
solvation of PE^ . 

Caleulated diffusion eoeffieients quantify the ionie motion in eaeh eleetrolyte and relate to 
the solvation struetures. We find that the largest Ei+ diffusion eoeffieient oeeurs in EMC. This 
is eonsistent with the measured viseosities of both organie solvents, and the diffusivity values 
obtained agree well with experimental values. The magnitude of the diffusion eoeffieient is largely 
influeneed by how tightly the ion is solvated, for both Ei+ and PE^ . A more tightly bound solvation 
strueture, sueh as for Ei+ in EC, leads to slower diffusion of the solvated ion. Comparison of Ei+ 
and PEg diffusion eoeffieients indieate that PEg diffuses faster than Ei+ in all the eleetrolytes 
examined here, even though it is the heavier speeies. This ean be attributed to the faet that the 
solvent moleeules internet more strongly with Ei+ than with PEg . 

Eurthermore, we quantified finite size and time seale effeets on the solvation struetures and 
diffusion eoeffieients using ReaxFF. We find that solvation struetures do not ehange signifieantly 
for larger system sizes and longer time seales than used here with DET, but there are some strue- 
tural differenees between ReaxFF and DET simulations. Absolute diffusion eoeffieients are more 
affeeted by size and time seale effeets, as uneertainties ean be as large as 50%, but relative values 
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remain consistent. 


Our work has shown that a more tightly bound solvation structure leads to slower diffusion, 
and a weakly bound solvation structure leads to faster diffusion. To improve the mobility of Li 
ions in solution, our results suggest that Li+ must have weak interactions with the organic solvent 
used in the electrolyte. This is valuable insight that can be used to improve the cycling rate of 
Li-ion batteries and potentially lead to the design of new electrolytes for better overall battery 
performance. 
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